/*  sdfrqrrt.cpp  by K.Tsuru  */
// function ID = 3013  DRADIX only
#include "sn.h"
static const char* const func = "Qrrt/RecQrrt";
/*****************************************************************
SDouble class
reciprocal quartic root of x
[Algorithm]
Get a solusion of an equation 1/(x^4) - a = 0 by the
Newton's method with "doubling effective figures"(DEF).
See Sqrt() for detail.
Let d(n) = x(n){1-a*x(n)^4}/4... faster than {x(n)-a*x(n)^5}/4 
because 1-a*x(n)^4 is very small
    x(n+1) = x(n) + d(n)
******************************************************************/
SDouble RecQrrt(const SDouble& a){
    if(a.Type() != a.REAL) a.SetError(a.RADIX_ERR, func, 3013);
    if(a.Sign() <= 0) a.SetError(a.DOMAIN_ERR, func, 3013);
    RealSize C;
    uint max_sz = a.MaxSize();
    SDouble w(Dabs(a)), one(1.0);
    double w0;

    if(a.IsOne()) return one;   // a = 1.0

    int e = w.NetRdxExp(), q, r, pow10 = 0;
// a = d*R^e = (d*R^q)*R^(4*r) , e = q + 4r
    r = e/4;
    q = 4*r;
    if(e - q > 0){ r++; q += 4; }
    w.MultPowRdx(-q);   // Here 1/R^4 <= w < 1.0.
/*set an initial value y0*/
    w0 = doubleD(w);
    while(w0 < 1.0e-4){ // 0.0001 <= w0 < 1.0
        w0 *= 10000.0;  pow10++;
    }
    if(pow10) w.MultPow10(4*pow10);
//enter the irreducible size mode
    SDouble y(a.REAL, max_sz), delta(a.REAL, max_sz), z(a.REAL, max_sz);
    y = pow(w0, -1.0/4.0);  // y has about DOUBLE_FIG figures.
/***************************************************************/
    int c = 0, itrmax = howpow2( (DFIGURES*max_sz)/DOUBLE_FIG+1 ) +6;
    uint ef = (DOUBLE_FIG*2u)/DFIGURES, fig = C.EffFigures();
    int fullPrec = 0;
    delta = w0 - w;
    if( delta.Sign() ){
    // If w0 is very close to double value, higher precision is necessary.
        ef = max( ef, 2u*(uint)abs( delta.NetRdxExp()) );
    }
//enter fiexd point mode
    y.FixedPoint(y.RdxExp());
    if(ef > fig) ef = fig;
    do{
        if( (ef = C.SetEffFig(ef)) >= fig) fullPrec = 1;
        z = y*y;    z = z*z;         // z = y^4
        delta = y*(one - w*z);
        delta = DsDiv(delta, 4);    // delta /= 4
        y += delta;
        c++;
        ef *= 2;
        if(ef >= fig) ef = fig;
        C.SetEffFig(0);
    } while( ( !delta.IsHidden() && (c < itrmax) ) || !fullPrec );
    y.iterationCount = c;
//#ifndef NDEBUG
    if(c == itrmax) y.SetError(y.FATAL, func, 3012);
//#endif
    y.PointFree();
    y.Reform(3011);

    if(y.Verify()){
        z = y*y;    z = z*z;    // z =y^4
        delta = one - z*w;  // delta = 1 - w*y^4
        if(!delta.IsMLT(w)){
            y.SetError(y.VERIFY, func, 3010);
        }
    }
    if(r) y.MultPowRdx(-r);
    if(pow10) y.MultPow10(pow10);
    return y;
}